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Abstract A low-order mimetic finite difference (MFD) method for Reissner-Mindlin 
plate problems is considered. Together with the source problem, the free vibration 
and the buckling problems are investigated. Full details about the scheme implemen- 
tation are provided, and the numerical results on several different types of meshes are 
reported. 



1 Introduction 

The Reissner-Mindlin theory is widely used to describe the bending behavior of an 
elastic plate loaded by a transverse force. However, its discretization by means of 
Galerkin methods is typically not straightforward. For instance, standard low-order 
finite element schemes exhibit a severe lack of convergence whenever the thickness 
is too small with respect to the other characteristic dimensions of the plate. This 
undesirable phenomenon, known as shear locking, is nowadays well understood: as 
the plate thickness tends to zero, the Reissner-Mindlin model enforces the Kirchhoff 
constraint, which is typically too severe for the discrete scheme at hand, especially if 
low-order polynomials are employed (see, for instance, the monograph by Brezzi and 
Fortin [16]). The root of the shear locking phenomenon is that the space of discrete 
functions which satisfy the Kirchhoff constraint is very small, and does not properly 
approximate a generic plate solution. The most popular way to overcome the shear 
locking phenomenon in Galerkin methods is to reduce the influence of the shear energy 

Lourenco Beirao da Veiga 

Dipartimento di Matematica "F. Enriques", Universita degli Studi di Milano, 
Via Saldini 50, 20133 Milano, Italy 
E-mail: lourenco.beirao@unimi.it 

Carlo Lovadina 

Dipartimento di Matematica, Universita di Pavia, 
Via Ferrata 1, 1-27100 Pavia, Italy 
E-mail: carlo.lovadina@unipv.it 

David Mora 

Departamento de Matematica, Universidad del Bi'o-BIo, Casilla 5-C, Concepcion, Chile and 
CI 2 MA, Universidad de Concepcion, Concepcion, Chile 
E-mail: dmora@ubiobio.cl 



2 



Lourengo Beirao da Veiga et al. 



by considering a (selective) reduced integration of the shear part, by resorting to a 
mixed formulation or by introducing a suitable shear reduction operator. Indeed, several 
families of methods have been rigorously shown to be free from locking and optimally 
convergent; let us mention, for instance, [1-4,14,17,28,29], and [37,38,40,43,44]. 

In the last years, many mimetic discretizations have been developed for the dis- 
cretization of problems in partial differential equations. The mimetic finite difference 
(or MFD) method has been successfully employed for solving problems of electromag- 
netism [35], gas dinamics [21], linear diffusion (see, e.g., [8,9,12,13,15,18,19,31,32,36] 
and the references therein), convection-diffusion [23], Stokes flow [6] and elasticity [5]. 
We also mention the development of a posteriori estimators for linear diffusion in [10] 
and post-processing technique in [22]. Finally, the mimetic discretization method has 
been shown to share strong similarities with the finite volume method in [26]. 

Recently, a mimetic finite difference (MFD) procedure has been proposed and the- 
oretically analysed for Reissner-Mindlin plates in [11]. The method, which can be con- 
sidered as a MFD version of the MITC and Duran-Liberman elements, combines the 
excellent convergence behaviour of the latter schemes with the great flexibility in han- 
dling the mesh of the former approach. The aim of this paper is to numerically assess 
the actual performance of the MFD method, by considering the source problem, as 
well as the free vibration and the plate buckling problems. 

A brief outline of the paper is as follows. In Section 2 we recall the Reissner- 
Mindlin plate problems, together with the necessary notations. Section 3 concerns with 
a presentation of the numerical scheme proposed and analysed in [11]. In particular, 
full details about the method implementation, that where missing in [11], are provided. 
In Section 4 we report the numerical results obtained using several types of meshes. 
We end the paper with some concluding remarks. 



2 The Reissner-Mindlin plate equations 

Here and thereafter we use the following operator notation for any tensor field r = (t^) 
i,j = 1,2, any vector field r) = (%) i = 1, 2 and any scalar field v: 

divrj := dir/i + d 2 rj2, rotrj := dir/ 2 - <9 2 ??i, tr(r) := m + r 2 2, 
y w ._ ( 9 ^ v \ curlu ■= ( ^ 2V \ divr ■= ( dlTl1 + d2l ~ 12 \ 

\d2Vj' ' \-dlVj' ' \dlT21 + d 2 T22j ' 

Consider an elastic plate of thickness t such that < t < diam(J7), with reference 
configuration SI x ( — 5, 5) , where SI is a convex polygonal domain of M 2 occupied 
by the midsection of the plate. The deformation of the plate is described by means 
of the Reissner-Mindlin model in terms of the rotations f3 = (/3i , /3a) of the fibers 
initially normal to the plate's midsurface, the scaled shear stresses 7 = (71, 72), and the 
transverse displacement w. Assuming that the plate is clamped on its whole boundary 
dSl, the following strong equations describe the plate's response to conveniently scaled 
transversal load g £ L 2 (S7): find (f3,w,f) such that 

' -divCe(/3) - 7 = in f2, 
-divf = g in SI, 

7 = Kt~ 2 (Vw-(3) in SI, [ > 

J 3 = 0,w = on9ft, 
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where the tensor of bending moduli is given by: 

CT := 12(1 E „2) ((1 " V)T + ^ tr(T)I) ' 

with E > representing the Young modulus, < v < 1/2 being the Poisson ratio for 
the material and I indicating the second order identity tensor. 
Let the Hq(Q) -elliptic bilinear form be given by 

a(/3,r?):= f Ce(/3) : e( V ) = E f [(1 - u)sifi) : s ( v ) + u div^divr?] , (2) 

with e = (sij)i<ij<2 the standard strain tensor defined by £ij(f3) ■= ^(dif3j + 
d j p i ),l<i,j<27' 

Then, the variational formulation of problem (1) reads: 

Problem 1 Find (/3,w) G Hl{of x Hq(O) such that 

a(A v) + Kt~ 2 (Vw - /3, V« - i|)o,n = (5, «) ,fl V(i|, w) G^ffifx tfp 1 («). 

In this expression, k := Efc/2(1 + 1/) is the shear modulus with k a correction factor 
usually taken as 5/6 for clamped plates. 

We will also consider the free vibration and buckling problem for plates. 

The free vibration problem of a plate is (see [25,27,28,30]): 

Problem 2 Find A G E and ^ (/3, w) G H^{nf x H%(f2) such that 
a(J3, 77) + nr 2 {Vw -0,Vv- 77)0/2 = A 

for all (tj, v) € ffo(J?) 2 x .Hq (,!?), where A = puj 2 /t 2 , with p being the density and w 
the angular vibration frequency of the plate and the corresponding eigenfunctions are 
the vibration modes. 

The buckling problem of a plate is (see [33,39]): 

Problem 3 Find \ bp G K and ^ {f3,w) G Hl(nf x Hl(Q) such that 

o(/3, i?) + Arf" 2 (Vw - /3, Vt) - ri) , n = \ bp (<rVw, Vv) ,a V(t|, v) G Hh {nf X H^(Q), 

where cr(x,y) G R 2x2 is a symmetric tensor which corresponds to a pre-existing stress 
state in the plate, \ bp = \ bc /t 2 , with X bc being the buckling coefficients of the plate 
and the corresponding eigenfunctions are the buckling modes. 

Accordingly with (1), for the problems above the scaled shear stresses can be com- 
puted by 7 = K,t~ 2 (Vw — f3). 

3 A Mimetic Finite Difference (MFD) discretization 

In this section we review the mimetic discretization method for the Reissner-Mindlin 
plate bending problem presented in [11], and extend it to the free vibration and buckling 
problems. Finally, in Section 3.6 we give the details on the implementation of the 
method. 
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3.1 Notation and assumptions 

Let {7~h}h be a sequence of decompositions of the computational domain S7 into N{Th) 
polygons E. We assume that this partition is conformal, i.e. intersection of two different 
elements E\ and E2 is either a few mesh vertices, or a few mesh edges (two adjacent 
elements may share more than one edge) or empty. We allow 7~h to contain non-convex 
and degenerate elements. For each polygon E, \E\ denotes its area, He denotes its 
diameter and h := max. Ee j- h ^E- 

We denote the set of mesh vertices and edges by Vh and Eh, the set of internal 
vertices and edges by V® and the set of vertices and edges of a particular element 
E by V/f and £^ , and the set of boundary vertices and edges by and £®, respectively. 
Moreover, we denote a generic mesh vertex by v, a generic edge by e and its length 
both by h e and |e|. 

A fixed orientation is also set for the mesh 7^, which is reflected by a unit normal 
vector n e , e £ £h, fixed once for all. Moreover, t e denotes the tangent vector defined 
as the counterclockwise rotation of n e by tt/2. 

For every polygon E and edge e £ £/f, we define a unit normal vector n|; that 
points outside of E, and by t E the tangent vector as the counterclockwise rotation of 
n e E by n/2. 

The mesh is assumed to satisfy the shape regularity properties detailed in [11]. 

We make also the following assumption on the data. The scalar functions E, v are 
piecewise constant with respect to the mesh 7/j. Moreover, there exist two positive 
constants C* and C* such that C* < E < C* on the whole domain. The above 
uniformity condition on E is standard, while the piecewise constant condition can be 
interpreted as an approximation of the data and is introduced only for simplicity. In 
the general case, it is sufficient to assume that E and v are (piecewise) W 1 ' 00 and to 
introduce an element-wise averaging in the data of the numerical scheme. 

3.2 Degrees of freedom and interpolation operators 

The discretization of Problems 1-3 requires to discretize the scalar field of displacement 
and the vector fields of rotations and shears. In order to do so, we introduce the degrees 
of freedom for the numerical solution in accordance with the correspondence 

w,v G Hq(O) -> w h ,v h 6 W h , 
0,rieHZ(nf ^/3 h ,r, h eH h , 

~f,5eL 2 (nf -^~f h ,5 h er h , 

where represents the linear space of discrete displacement, indicates the linear 
space of discrete rotations and r h is the linear space of discrete shears. 

The discrete space for transverse displacements is defined as follows: a vector 
u/j G Wh consists of a collection of degrees of freedom 

v h ■= {v V } veV ° h -< 

one per internal mesh vertex, e.g. to every vertex v £ VjJ, we associate a real number 
v v . The scalar v v represents the nodal value of the underlying discrete scalar field of 
displacement. The number of unknowns is equal to the number of internal vertices. 
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The discrete space for rotations is defined as follows: a vector r) h £ Hy t is a 
collection of degrees of freedom 

Vh = {y}vev°> 

i.e. we assign a vector if £ M 2 per each vertex v £ V°. The vector rf represents 
the nodal values of the underlying discrete vector field of rotations. The number of 
unknowns is equal to twice the number of internal vertices. 

Finally, the space for the discrete shear force i"), is defined as follows: to every 
element E in 7^ and every edge e £ £/f PI £ °, we associate a number S%, i.e. 

We make the continuity assumption that for each edge e shared by two element E\ 
and Ei, we have 

3%! = Se 2 - 

The scalar 8^ represents the average on edges of the discrete shears in the tangential 
direction. The number of unknowns is equal to the number of internal edges. 

We now define the following interpolation operators from the spaces of smooth 
enough functions to the discrete spaces Wh, and F^, respectively. For every function 
v £ C°(n) n H^(O), we define v t £ W h by 

v{:=v(v) Vv£V£. (3) 
For every function r) £ [C°(i?) H H^(f2)] 2 , we define r/i £ by 

Vi--=t,{m) Vv£V,°. (4) 
For every function S £ i?o(rot; Q) n L s (f}) 2 , s > 2, we define <5n £ /ft by 

(S u )e ■= ^ lj ■ t% VEe% Ve££fn£t (5) 

For all E £ 7/j in the sequel we will also make use of local interpolation opera- 
tors i>i e, Vi e> &HE> with values in W^Ie, H^e, E^e respectively; such operators are 
simply the obvious restriction of the global ones to the element E for functions which 
are sufficiently regular on E. 

Remark 1 We note that in the present paper we are considering the scheme of [11] 
without the edge bubbles, see Remark 4 of [11]. Such version of the method is more 
efficient in terms of accuracy vs number of degrees of freedom, while the loss of stability 
is seen only on very particular mesh patterns. Indeed, in the numerical test of Section 4, 
only the first family of (triangular) meshes suffers from such drawback. 
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3.3 Discrete norms and operators 
We endow the space with the following norm 

IKIlk := E W v h\\w h ,E = E I s ! E 

EeT h EeT h eeef 

where vi and V2 are the vertices of e. Although irrelevant in (6), in the following we 
will always consider that vi and V2, the vertices of a generic edge e, are oriented in 
such a way that t e E points from vi to V2 . 
In the space H^, we consider the norm 

\\\v h \\\ 2 Hh ■= E \\\r>h\\\ 2 H h ,E= E i^i E 

-EeTh EeT h eeff 

where vi and V2 are the vertices of the edge e, and 1 1 ■ 1 1 denotes the euclidean norm on 
vectors. 

In the space 7\, we consider the following norm 

\\s h \\r h ■■= E w s h\\r h ,E= E i £ i E i^i 2 - ( g ) 

EeTh EeT h e££Z 

The norms on and are H 1 (f2) type discrete semi-norms, which become norms 
due to the boundary conditions on the spaces, while the norm for i^i is an L 2 (f2) type 
discrete norm. 

In the sequel we will also use the following norm on H^, which is a ||e(-)||o Q type 
discrete norm: 

\\Vh\\H h ■= E \\nh\\ 2 H h ,E= E ™|||ry h - c([-y,s]) Ij£; |||l fhijB , (9) 
EeTh EeTh C 

where (x, y) are local cartesian coordinates on E which are null on the barycenter of E, 
so that the function [— y, x] represents a (linearized) rotation around the barycenter. 

We now introduce the discrete gradient operator V/j, defined from the set of nodal 
unknowns Wh to the set of edge unknowns r), as follows: 

V,, : w h -> r h 

where vi and V2 are the vertices of e. 

We consider also a reduction operator, defined from the discrete space of rotations 
Hh to the set of edge unknowns i"), as follows: 

n h ■. H h -> r h 

(i7 ft T7 h )|:=i[r, Vl +77 V2 ]-t| VEeT h , Ve g £jf , Vi| h g H h , 
where vi and V2 are the vertices of e. 



V 1 ) 



(6) 



(7) 
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3.4 Scalar products and bilinear forms 

We equip the space F^ with a suitable scalar product, defined as follows: 

hh> S h]r h ■= b/h> S h]r h ,E, (10) 
EeT h 

where [■, -]r h .E is a discrete scalar product on the element E. 

The scalar product must satisfy the following stability and consistency conditions 
(see [11]). 

(SI) There exist two positive constants c\ and ci independent of h such that, for 
every 5^ £ F^ and each E £ 7/j, we have 

ci||*/ 1 ||j> 1 ,e < [Sh,Sh\r h ,E < c2\\d h \\r htE . (11) 



(S2) For every element E, every scalar linear function pi on E and every 5^ 6 F/j, 
we have 

[(curlpi) u ,6 h ]r h ,E = / Pi(rotr h 8 h ) E - V 5% / pi, (12) 

J E ^ Je 

where the operator (rot^ (i S^)^ := jj^ X^ee£ E ^sl e l' 

We denote with a^(-, ■) : x — > R the discretization of the bilinear form 
et(-, ■), defined as follows (see (2)): 

ah(Ph, , nh)= Y ah(Ph,Vh) VPh,Vh€ H h, (13) 

EeT h 

where af'(-, •) is a symmetric bilinear form on each element E, mimicking 
a h{PinVh) ~ / Ce(/3 fc ) : e(ij h ) . 



JE 

Similarly to the previous case, we introduce a stability and consistency assumptions 
for the local bilinear form af^(-, •). 

(Sl a ) there exist two positive constants ci and C2 independent of h such that, for 
every %£ and each E £ Th, we have 

ci||»7hllff h ,.E! < a h(Vh,Vh) < C2\\Vh\\H h ,E- (14) 

(S2 a ) For every element E, every linear vector function pi on E, and every r] h £ H^, 
it holds 



E 



((pi)l,i?fc) = E [(Ce(pi)n^-(f I^W 2 ])]- (15) 



Remark 2 The scalar product and the bilinear form shown in this section can be built 
element by element in a simple algebraic way. The details are shown in Section 3.6. 
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3.5 The discrete method 

Finally, we are able to define the mimetic discrete method for Reissner-Mindlin plates 
proposed in [11]. Let the loading term 

{g,v h ) h := ^ 9\e^v V 'u e , (16) 
E£T h i=l 



where vi , . . . , v mE are the vertices of E, g\ E '■= pjy J E g, and oj e , . . . , u E XB are positive 
weights such that 

rriE 

Pi = y^Pl(vj)^g 



E i=l 



for all linear functions p\ . The loading term above is an approximation of 

{g,v h ) h ~ / gv. 



J a 

Then, the discretization of Problem 1 reads: 
Method 1 Given g £ L 2 {Q), find (f3 h ,w h ) £ H h x Wh such that 

o-hiPh^h) + nt~ 2 [V h w h - II h (3 h ,X7 h v h - n h -n h ]r h = (g,v h ) h 

for all {r} h ,v h ) G H h x W h . 

In order to extend the method to the free vibration problem, we introduce the 
following mass bilinear form in Hh x Wh 

m h (l3 h ,w h ;r] h ,v h ) = ^ mf l {p h ,w h \r) h ,v h ) (17) 
EeTh 

for all Phy'Hh £ Hh and u>h,Vh £ Wh, where the local forms 

rriE ^2 m E 

™h(0h> w h;vh> v h) = ^2 wVivV '^E + j2 X^ Vl '^'We- 

1=1 1=1 
Then, the discretization of Problem 2 reads: 
Method 2 Find A/j £ R and {f3 h iWh) £ Hh x Wh such that 

o-hiPh^h) + Kt~ 2 [V h w h - II h (3 h ,V h v h - n h r) h \r h = A/j m h (/3 h ,w h ;ri h ,v h ) 

for all (r]h,vh) £ H h x W h . 

Finally, in order to discretize the buckling problem we introduce a discrete bilinear 
form 

bh{w h ,v h ) = ^2 bh(w h ,v h ) Vu>h,v h eWh, (18) 
EeT h 

where ■) is a symmetric bilinear form on each element E, mimicking 

bh(wh,v h ) ~ / (o-Vw h ) ■ Vv h ■ 

JE 

We assume for simplicity that the stress datum <j is piecewise constant on the mesh, a 
condition that can also be considered as an approximation of a given data. We require 
that the local bilinear forms 6/f(-, •) satisfy the following stability and consistency 
conditions. 
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(Slf,) There exists a positive constant c independent of h such that, for every G Wh 
and each E G Th, we have 

bf (vh, v h) < c\\v h \\w h ,E- ( 19 ) 

(S2&) For every element E, every scalar linear function pi on E, and every G Wh, 
it holds 

6f ((pi)i,« h ) = £ (<rVpi ■n|)M[« Vl + ,/*], (20) 

where we recall that cr\ E G R 2x2 is constant and symmetric. 

Such a condition asserts that the discrete bilinear form is exact when tested on linear 
functions. We also remark that for the form (•, ■) we do not require any lower bound, 
such as the one in (14). Indeed, assuming a lower bound condition for b^ (•, •) would 
be unnatural, since the stress datum a can be a singular second-order tensor. 
The discretization of Problem 3 then reads: 

Method 3 Find X h h p G R and (f3 h ,w h ) G H h x W h such that 

a-hiPh, Vh) + Kt~ 2 [^ h w h - II h (3 h , V h v h - n h T) h ]r h = b h (w h , v h ) 
for all (f] h ,v h ) G H h x W h . 



3.6 Implementation of the method 

In this section we describe explicitly how to build the local bilinear forms appearing in 
the previous sections. 

In what follows m = m(E) G N will indicate the number of vertices of the polygon 
E. We number the vertices in counterclockwise sense as vi, .., v m and analogously for 
the edg es ei , .. , em , so that \jj and Vj+i are the endpoints of edge ej, j = 1,2, ..,m. 
Note that here and in the sequel all such indexes are considered modulus m, so that 
the index m + 1 is identified with the index 1. There are a total of 3m local degrees of 
freedom associated to each element of the mesh, three for each vertex. We order such 
local degrees of freedom first with all rotations and then all deflections, ordered as the 
vertices 

{Ve'Ve'-'^E , v E> v E>-< v E h 

where {rj E ,v E ) G H h \ E x W h \ E - 

The final local bilinear forms M = M(E) G R 3mx3m associated to each element E 
will be the sum of two parts 

M = Mi + ^~ 2 M 2 , (21) 

the first one being associated to the a^(-,-) term and the second one to the shear 
energy term. Once the elemental matrices M are built, the global stiffness matrix is 
implemented with a standard assembly procedure as in classical finite elements. 
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3.6.1 Matrix for the bilinear form a^(-, •) 

We start from the bilinear form a^(-, •), which is the sum of local bilinear forms that 
we express as matrices M = M(E) G R 2mx2m 

af l {(3 E ,r lE )=[3 T E Wr lE \/E e T h , V0 E , Ve € H h \ E . 

The first and main step is to build the matrix M. With this purpose we introduce the 
matrices N = N(E) and R = R(E) in R 2mx6 . Note again that for ease of notation we 
do not make explicit the dependence on the involved matrices from E. Let qi, . . . , q@ be 
the following basis for the first order vector polynomials (with 2 components) defined 
on E: 




Then, the six columns IMi, . . . , Ne of N are vectors in R 2m defined by the interpolation 
of the polynomials qi, . . . ,qg into the space Hf l \ E (see (4)) 

Nj = (q.j)i,E 

so that for 1 < i < m and 1 < j < 6 

The columns of the matrix N thus represent the linear polynomials q^ written in terms 
of the degrees of freedom of H] l \ E . 

The columns Rj of the matrix R are instead defined as the vectors in R associated 
to the right hand side of the consistency condition (S2 a ), computed with respect to 
the polynomials qj, j = 1, . . . , 6. In other words Rj is the unique vector in M 2m such 
that for all r} E £H h \ E = R 2m 

771 . . 

(RjfvE = E ( Ce (<i> e i) ■ (^Mg + vT 1 }) 

i=l 

see equation (S2 a ). Note that, since e(qj) = for j = 1,2,3, the first three columns 
Ri, R2, R3 of R have all zero entries. 

From the definition of the vectors and Rj, it is clear that the consistency con- 
dition (S2 a ) translates into the algebraic condition 

MNj = Rj j = l,...,6 <^ MN = R. (22) 

We therefore introduce the matrix K £ R 6x6 defined by 

K = N T R = R T N. 

It is easy to check that such matrix is symmetric and semi-positive definite. Moreover, 
it is of the form 

K = I 03x3 ^3x3^ 
~ V° 3X3 K * / 



Numerical results for mimetic discretization of Reissner-Mindlin plate problems 



11 



with K* positive definite. Therefore is it immediate to compute the pseudo inverse of 
K 



K 1 



t _ / 03x3 03x3 

3 x3 K- 1 



We are now ready to define the local matrix M. Let P be a projection on the space 
orthogonal to the columns of N 

P = I 2 mx 2m -N(N T N)- 1 N T 

with I2 m x2m the identity matrix. We then set 

M = RK f R T + aP 

with a 6 R any positive number, typically scaled as the trace of the first part of the 
matrix. Then, it is immediate to check that M satisfies the consistency condition (22). 
Moreover, the positivity up to the kernel is easy to check, while the uniform positivity 
represented by the stability condition (Sl a ) can be proved with the techniques shown 
in [20,6]. 

Finally, note that the matrix M £ R 2mx2?T1 i s defined only with respect to the rota- 
tion degrees of freedom, since the bilinear form a^(-, •) is independent of the deflection 
variable. When it comes to build the local matrix Mi £ M 3 " 7,x m appearing in (21) 
one simply needs to introduce the restriction matrix S € ^ 3mx2m 



"rax2m 



SMS T . 



S 

and set 

I 

3.6.2 Matrix for the shear term 

The local matrices for the shear part are obtained as a product of matrices representing 
the operators and bilinear forms that appear in the second term of the left hand side 
of Method 1. We therefore start building a matrix M~ = M(E) £ E mXm that represents 
the local scalar product 

bfE> s E]r H ,E = nrEM6 B v~i E ,s E e r h \ E . 

We order the m degrees of freedom of E^\ E as the edges of E. The construction follows 
the same philosophy as in the previous section and therefore is presented more briefly. 
Now, the two columns of the matrix N £ R mx2 are defined by 

Hj = (curl qj)n >E j = 1,2, 

where the sub- index n represents the interpolation operator shown in (5) and qi , qi 
denote the following basis of the (zero average) linear polynomials on E 

qi =x, q 2 = y. 

Analogously, the matrix R £ R mx2 is defined through its columns as the right hand 
side of (S2) 

(Rj) T s E = -Y, s e / m Vj = 1, 2, vs E £ r h \ E = u m , 
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where we neglected the rotr h part since qi and q 2 have zero average on E. Again, we 



need to introduce K G 



given by 



K = N T R = R T N 



R(K)~ 1 R T + aP 



that is easily shown to be positive definite and symmetric. We can therefore finally set 

M 

with q G M + and the projection matrix 

P = Imxm " 



1 xm-N(N T N)- 1 N T . 

The consistency condition MN = R follows by construction while the stability can be 
derived with the results in [20]. 

The local matrix M2 appearing in (21) can be built combining M with a matrix 
C = C(E) G jjraxsm re p regen t} n g y h anc j jj h operators that appear in Method 1. 
We therefore set 

C= (-Ci C a ) 
with the matrix Ci = d(E) G K mx2m 



representing the II ^ operator 



/(**;) (tljr «ix2 1X2 lx2 \ 

0ix2 (tg) T (tg) T 1X2 1X2 



01x2 o lx2 (t e j) J {t^y Oi 



\(t^) T 0i x2 oi x2 oi x2 (t^)v 

and the matrix C 2 = C 2 (75) G R mxm representing the V/j operator 





1- 


leiP 1 


leiP 1 




















-leap 1 


leap 1 










c 2 = 










-leal" 1 


\e3\- 1 









{- 


6j7t 1 










|e,„r 


-V 



Finally, the local matrices for the shear part are given by 



C T MC. 



3.6.3 Right hand sides 

The loading term for the source problem in Method 1 follows immediately from (16). 
One gets the local right hand vectors b = b(E) G R 3m defined by 







if j = 1,2, ..,2m 
if j = 2m + 1,2m + 2,. .,3m, 

that are then assembled as usual into the global load vector. 



-1 (j— 2m) 
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The mass matrix for the free vibration problem in Method 2, associated to the 
bilinear form (17) is built again by a standard assembly procedure. The local mass 
matrices D = D(E) G R 3mx3m associated to the elemental mass bilinear forms 

m h(PE' w E;VE> v E) = (Pe> w e) t d {Ve^e) V£ G %, 



VPe'Ve € H h \ E , \/w Ey v E G W h \ E 



are diagonal and defined by 



t 2 w^ /2l /12 if i= 1,2, ..,2m 

(i-2m) 



Hi = 2m + 1,2m + 2, ..,3m 



where the symbol \ ] stands for a round up to the nearest integer. 

The stress matrix for the buckling problem in Method 3, associated to the bilinear 
form (18) is also built as the sum of local matrices B = B(E) G R mXm 

bh ( w £.»e) = w E B v E \/w E ,v E G W h \ E . 

Note that the symmetric tensor er G KL 2x2 that appears in (S2j) can have either rank 
2 or rank 1. In order to build the matrix B, we start introducing {51,(72,93} a basis 
for the linear polynomials on E, such that q\ = 1 and 52, <?3 have zero integral on E. 
Moreover, if rank(cr) = 1, we also require that V<j2 G Ker(tr). We then define as usual 
the auxiliary matrices N = H{E) G E mx3 and R = R(E) G R mx3 through its columns. 
We set 

N j = (q j )i,E, 3 = 1,2,3, 

where the sub- index 1 denotes the interpolation operator in (3), and define Rj as the 
unique vector in K m such that 

rn . 

*jVE = Y, ( ctV ^ ■ n S) ^'E + v e +1 \ V i = 1. 2- 3- 6 W h \ E = WL m 



in accordance with (S2{,). Note that clearly Ri is null, and that, if rank(tr) = 1 also 
R2 is null. One then defines as usual the semi-positive definite and symmetric matrix 
R = K(E) G K 3x3 

R = R T N = N T R. 

Since K is block diagonal, with the first block of zeros and the second invertible, it is 
immediate to compute the pseudo inverse matrix , in a way similar to the one used 
for K in Section 3.6.1. Then, we introduce B = B(E) G K mxm 

B = R (K) f R T + aP 

with a£l non negative and the projection matrix P = I mxm — N(N T N) _1 N T . Note 
that, since no global coercivity conditions are required, differently from the previous 
matrices also the choice a = can be taken. 

Finally, note that the matrix B G R mxm is defined only with respect to the deflec- 
tion degrees of freedom, since the bilinear form bf l (-,-) is independent of the rotation 
variable. The remaining entries in the assembled (right hand side) stress matrix asso- 
ciated to Method 3 can be simply filled with zeros. 
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4 Numerical results 

The numerical method analyzed has been implemented in a MATLAB code. 

For all the computations we took Q : = (0, l) 2 , for the Young modulus we choose: 
£7 = 1. 

We have tested the method by using different meshes. We report the results ob- 
tained using the families of meshes shown in Figure 1 to Figure 7. 



— 71 : Triangular mesh. 

— T h : Structured hexagonal meshes. 

— 7/f : Non-structured hexagonal meshes made of convex hexagons. 

— 7/f : Regular subdivisions of the domain in N x N subsquares. 

— 7/f: Trapezoidal meshes which consist of partitions of the domain into N x N 
congruent trapezoids, all similar to the trapezoid with vertices (0, 0), 0), (^, |), 
and (0, |). 

— 7/f : Regular polygonal meshes built from considering the middle point of each 
edge as a new node on the mesh; note that each element has 6 edges. 

— 7^: Irregular polygonal meshes built from 7/f moving randomly the middle point 
of each edge; note that these meshes contain non-convex elements. 




N = 4 

Fig. 1 Square plate: meshes 7?. 



N = 8 
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N = 4 

Fig. 4 Square plate: meshes fc- 



N = 8 
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N = 4 

Fig. 7 Square plate: meshes . 



N = 8 
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4.1 Source problem 

As a test problem we have taken an isotropic and homogeneous plate, clamped on the 
whole boundary, for which the analytical solution is explicitly known (see [24]). 
Choosing the transversal load g as: 

g(x, y) = 12{1 E _ U 2) [ 12 y(y - ^ 5x2 - 5x + x ) O 2 ^ - x ) 2 + x( - x - v^y 2 ~ 5 y + 

+12x(x - l)(5y 2 - 5y + 1) (2x 2 (x - l) 2 + y(y - l)(5x 2 - 5x + 1)) 1, 
the exact solution of the problem is given by: 



w(x,y) = ^x 3 (x - l) 3 y 3 {y ~ I) 3 

- — \y 3 {y - lfx{x - l)(5x 2 - 5x + 1) + x 3 {x - l) 3 y(y - l)(5y 2 -5y + l)\, 

5(1 — u) L J 

PAx.v) = v s (v - ifx 2 {x - \) 2 {2x - i), 

f3 2 (x,y) = x 3 (x - l) 3 y 2 (y - l) 2 (2y - 1). 

We have used a Poisson ratio v = 0.3 and a correction factor k = 5/6. 
The convergence rates for the transverse displacement w and rotations f3 axe shown 
in the following norms: 

e(pjo •— r-^H ' e \ w )o ■— i — i , (23J 

e (P h ■— 77; — „ m it, > e{w)i .— -jz . 

^(/3i,/3i) 1/2 [V^V^f 

(24) 

In (24) the bilinear forms a^(-,-) and [v]/ 1 ?, are exactly the ones denned in (13) 
and (10), respectively. We notice that it holds: 

\\Pi\\H h \\ w i\\w h 

Therefore, (23) and (24) represent discrete L°° and H 1 relative errors, respectively. 

Also, we define the experimental rates of convergence (rc) for the errors e(/3) and 
e(w) by 

log(e(-)/e'(-)) 



=(•) 



log(tyfc') 



where h and hi denote two consecutive meshsizes and e and e', respectively, denote 
the corresponding errors. 

Table 1 shows the convergence history of the Method 1 applied to our test problem 
with four different family of meshes. Table 2 shows instead an analysis for various 
thicknesses in order to assess the locking free nature of the proposed method. 
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Table 1 Convergence analysis for t = 0.01. Errors and experimental rates of convergence for 
variables (3 and w. 



Mesh 


l/h 


e(/3)o 


rc(/3) 


e(iu) 


rc(w) 


e(/3)i 


rc(/3)i 


e(iu)i 


rc(w) i 




8 


5.018o-2 




2.641c-2 




9.700o-2 




6.480c-2 






16 


1.103e-2 


2.19 


9.522e-3 


1.47 


2.967e-2 


1.71 


1.988e-2 


1.70 


T 2 

'h 


32 


2.788e-3 


1.98 


2.761e-3 


1.79 


8.992e-3 


1.72 


5.404e-3 


1.88 




64 


7.166e-4 


1.96 


7.351e-4 


1.91 


2.886e-3 


1.64 


1.403e-3 


1.95 




128 


1.814o-4 


1.98 


1.891c-4 


1.96 


l.OlOc-3 


1.51 


3.572c-4 


1.97 




8 


7.137o-2 




5.208c-2 




1.325o-l 




8.680c-2 






16 


3.505e-2 


1.03 


2.095e-2 


1.31 


5.465e-2 


1.28 


3.303e-2 


1.39 


T 3 

'h 


32 


1.1310-2 


1.63 


6.219e-3 


1.75 


1.793e-2 


1.61 


9.714e-3 


1.77 




64 


3.108e-3 


1.86 


1.634e-3 


1.93 


5.619e-3 


1.67 


2.571e-3 


1.92 




128 


7.991o-4 


1.96 


4.156c-4 


1.98 


1.846o-3 


1.61 


6.571c-4 


1.97 




8 


3.224o-2 




6.519c-2 




4.370o-2 




9.599o-2 






16 


8.156e-3 


1.98 


1.605e-2 


2.02 


1.132e-2 


1.95 


2.518e-2 


1.93 


T 4 

'h 


32 


2.051e-3 


1.99 


3.997e-3 


2.00 


2.866e-3 


1.98 


6.365e-3 


1.98 




64 


5.138e-4 


1.99 


9.983c-4 


2.00 


7.188e-4 


2.00 


1.595e-3 


2.00 




128 


1.285o-4 


1.99 


2.496c-4 


2.00 


1.798o-4 


2.00 


3.991o-4 


2.00 




8 


7.190o-2 




1.057c-l 




1.949o-l 




1.318c-l 






16 


1.677e-2 


2.10 


2.331e-2 


2.18 


1.127e-l 


0.79 


4.339e-2 


1.60 




32 


3.509e-3 


2.26 


5.201e-3 


2.16 


4.213e-2 


1.42 


1.080e-2 


2.01 


64 


5.942e-4 


2.56 


1.221e-3 


2.09 


1.127o-2 


1.90 


2.380e-3 


2.18 




128 


1.504o-4 


1.98 


2.896e-4 


2.08 


3.802o-3 


1.57 


5.516o-4 


2.11 



We observe from Table 1 that a clear rate of convergence 0(h 2 ) is attained for j3 and 
w for all family of meshes in the discrete L°° norm. Moreover, a rate of convergence 
0(h 3 / 2 ) for /3 and 0(h 2 ) for w for all family of meshes in the discrete H 1 norm 
have been obtained. Actually, the computation of (3 using meshes 7/f seems to be 
superconvergent . 



Table 2 Locking-free analysis for variable w (e(io)i) 



Mesh 


l/h 


i=1.0e-2 


t=1.0e-3 


t=1.0e-4 


t=1.0e-5 




8 


2.040179e-l 


9.381597e-l 


9.993307e-l 


9.999933e-l 


' h 


16 


2.975046e-2 


3.790016e-l 


9.773959e-l 


9.997674e-l 


32 


4.320781e-3 


4.271249e-2 


5.521358e-l 


9.902023e-l 




64 


8.636150e-4 


5.651503e-3 


8.549716e-2 


7.775880e-l 




8 


8.680034e-2 


8.648290e-2 


8.647974e-2 


8.647905e-2 


n 


16 


3.303805e-2 


3.289787e-2 


3.289649e-2 


3.289827e-2 




32 


9.714549e-3 


9.670538e-3 


9.670075e-3 


9.671654e-3 




64 


2.571361e-3 


2.558931e-3 


2.558797e-3 


2.555879e-3 




8 


9.598706e-2 


9.605054e-2 


9.605176e-2 


9.605118e-2 




16 


2.518265e-2 


2.518620e-2 


2.518623e-2 


2.518663e-2 




32 


6.364552e-3 


6.364126e-3 


6.364108e-3 


6.364647e-3 




64 


1.595350e-3 


1.595135e-3 


1.595115e-3 


1.595949e-3 




8 


7.471495e-2 


7.399270e-2 


7.398561e-2 


7.398560e-2 




16 


2.223311e-2 


2.217105e-2 


2.217068e-2 


2.217124e-2 


32 


6.155689e-3 


6.173926e-3 


6.174458e-3 


6.176650e-3 




64 


1.268306e-3 


1.308765e-3 


1.309856e-3 


1.311575e-3 




8 


8.776844e-2 


8.835775e-2 


8.547591e-2 


8.870029e-2 


n 


16 


3.128646e-2 


3.002581e-2 


3.088068e-2 


2.990004e-2 


32 


8.776518e-3 


8.804238e-3 


9.034351e-3 


8.770496e-3 




64 


1.925048e-3 


2.025191e-3 


2.042556e-3 


1.966998e-3 
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We observe from Table 2 that our Method 1 lead to wrong result only for triangular 
meshes 7^ , when the thickness of the plate is small. For any other family of meshes 
the method is locking-free. We also note that adding the middle point of each edge as 
a new node on any triangular mesh, the method is locking free, see row corresponding 
to Tl and 7jf. 

Remark 3 We note that the different behavior among the triangular mesh 7/f and the 
remaining grids is not surprising. Indeed, resembles a plain PI element, that is 
known to suffer from locking in the plate finite element literature, unless edge bubbles 
are added to the rotations (see also Remark 1). Moreover, the analysis of [11] does not 
apply to the 7/f meshes since the P1/P0 element is not stable for the Stokes problem. 
Instead, meshes Tff, 7jf fall into the hypotheses of the convergence theorem of [7], 
see Remark 4 of [11]. Meshes 7/j\ 7/f have a strong connection with the MITC4 finite 
element for plates that is known to be stable. Finally, meshes 7/f , 7^ again fall into 
the convergent cases considered in [7] and thus stable also for Reissner-Mindlin (see 



4.2 Free vibration of plates 

The effectiveness of the MDF method for free vibration analysis are demostrated by 
examples with different thickness and different boundary conditions. 

We have computed approximations of the free vibration angular frequencies u) h = 



here to mn are the computed frequencies, where m and n are the numbers of half- waves 
in the modal shapes in the x and y directions, respectively. L is the plate side length. 

We have considered a square plate of side length L = 1 and p = 1 and three 
different thickness t = 0.1, t = 0.01 and f=1.0e-5. We have also considered three 
different types of boundary conditions: a clamped plate (denote by CCCC), a simply 
supported plate (denote by SSSS), and a plate with a free edge (with three clamped 
edges and the fourth free, we denote by CCCF). 

In the following numerical tests, we show the results for the four lowest vibration 
frequencies. We tested also higher frequencies with similar results. 

Tables 3 and 4 show the four lowest vibration frequencies computed by Method 2 
with successively refined meshes of each type for a clamped plate with thickness t = 
0.1 and t = 0.01, respectively. The table includes orders of convegence, as well as 
accurate values extrapolated by means of a least-squares fitting. Furthermore, the last 
two columns show the results reported in [25,27,30]. In every case, we have used a 
Poisson ratio v = 0.3 and a correction factor k = 0.8601. The reported non-dimensional 
frequencies are independent of the remaining geometrical and physical parameters, 
except for the thickness-to-span ratio. 



[11])- 



ty ^k. In order to compare our results with those in [25,27,28,30], 
frequency parameter is defined as: 



a non-dimensional 
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Table 3 Lowest non-dimensional vibration frequencies for a CCCC square plate and t = 0.1 



Mesh 


Mode 


Af = 32 


TV = 64 


jV = 128 


i"dci' 


IjAIiI ClJJ . 


ho 


27 






1.5938 


1.5918 


1.5912 


1.84 


1.5910 


1.591 


1.5910 


T- 


U321 








1 .OO 


o.ujoy 


o.uoy 






w 12 


^ n^nn 

O. UtiUU 


3.0419 


3.0397 


1.90 


o.uooy 


3.039 


O.l/OOO 






A 9£07 






1 .OJ 


4 2624 




4 2624 






1.5958 


1.5923 


1.5914 


1.84 


1.5910 


1.591 


1.5910 


n 




3.0524 


3.0426 


3.0399 


1.83 


3.0388 


3.039 


3.0388 




^12 


3.0570 


3.0438 


3.0402 


1.88 


3.0388 


3.039 


3.0388 




<^22 


4.2903 


4.2701 


4.2645 


1.84 


4.2623 


4.263 


4.2624 




Wll 


1.5961 


1.5923 


1.5914 


1.97 


1.5910 


1.591 


1.5910 


T 4 

'h 


UJ21 


3.0526 


3.0424 


3.0398 


1.98 


3.0389 


3.039 


3.0388 




W12 


3.0526 


3.0424 


3.0398 


1.98 


3.0389 


3.039 


3.0388 




<^22 


4.2914 


4.2699 


4.2644 


1.97 


4.2625 


4.263 


4.2624 




Wll 


1.5967 


1.5925 


1.5914 


1.98 


1.5910 


1.591 


1.5910 


T 5 


^21 


3.0527 


3.0424 


3.0398 


1.98 


3.0389 


3.039 


3.0388 




t*>12 


3.0573 


3.0435 


3.0401 


2.00 


3.0389 


3.039 


3.0388 




1^22 


4.2943 


4.2705 


4.2645 


1.98 


4.2625 


4.263 


4.2624 



Table 4 Lowest non-dimensional vibration frequencies for a CCCC square plate and t = 0.01 



Mesh 


Mode 


AT = 32 


N = 64 


N = 128 


Order 


Extrap. 


[30] 


[25] 






0.1757 


0.1755 


0.1754 


1.89 


0.1754 


0.1754 


0.1754 




1^21 


0.3582 


0.3576 


0.3574 


1.87 


0.3574 


0.3574 


0.3576 




1^12 


0.3587 


0.3577 


0.3575 


1.91 


0.3574 


0.3574 


0.3576 




^22 


0.5289 


0.5272 


0.5267 


1.86 


0.5265 


0.5264 


0.5274 






0.1759 


0.1755 


0.1754 


1.87 


0.1754 


0.1754 


0.1754 


n 


^21 


0.3590 


0.3578 


0.3575 


1.84 


0.3574 


0.3574 


0.3576 




1^12 


0.3596 


0.3580 


0.3575 


1.87 


0.3574 


0.3574 


0.3576 




1^22 


0.5304 


0.5276 


0.5268 


1.83 


0.5265 


0.5264 


0.5274 




Wii 


0.1759 


0.1755 


0.1754 


1.99 


0.1754 


0.1754 


0.1754 




^21 


0.3593 


0.3579 


0.3575 


2.00 


0.3574 


0.3574 


0.3576 




^12 


0.3593 


0.3579 


0.3575 


2.00 


0.3574 


0.3574 


0.3576 




^22 


0.5306 


0.5275 


0.5268 


1.99 


0.5265 


0.5264 


0.5274 






0.1762 


0.1756 


0.1754 


2.21 


0.1754 


0.1754 


0.1754 


n 


1^21 


0.3597 


0.3579 


0.3575 


2.10 


0.3574 


0.3574 


0.3576 




^12 


0.3613 


0.3582 


0.3576 


2.33 


0.3574 


0.3574 


0.3576 




1^22 


0.5323 


0.5278 


0.5268 


2.16 


0.5265 


0.5264 


0.5274 



It can be seen from Tables 3 and 4 that our method converges with a quadratic 
order. 

Table 5 shows the four lowest vibration frequencies computed by Method 2 with 
successively refined meshes of each type for a clamped plate with £=1.0e-5. The table 
includes orders of convegence, as well as accurate values extrapolated by means of a 
least-squares fitting. In every case, we have used a Poisson ratio v = 0.3 and a correction 
factor k = 0.8601. The reported non-dimensional frequencies are independent of the 
remaining geometrical and physical parameters, except for the thickness-to-span ratio. 
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Table 5 Lowest non-dimensional vibration frequencies for a CCCC square plate and t=1.0e-5 



Mesh 


Mode 


N = 8 


N = 16 


N = 32 


Oi'der 


IjAIiI ClJJ . 




CJH 


fl 7653e-1 


0.1289e-l 


f) 1987e-2 


2.54 


-0 2995e-3 


-r 1 
>h 


Co' 21 


n 1 Q79^ n 
u. ly / ze-u 


u. ozzoe- r 


ft AQA9p 9 


9 


-u.joyuc-o 




^12 


n 99^np n 

17. iiOUc U 


zLflfifle 1 
u.^iuuue- r 


17. JUOOc-Z 


2.36 


fl 3^1 Qe 9 




^22 


n ^qi n 

u.oy ijc-u 


u.uuyoc-i 


u.yoyoL-z 


9 7fl 


u. ( uoue-o 




CJll 


0.1848e-3 


0.1778e-3 


0.1761e-3 


2.04 


0.1756e-3 




W21 


0.3927e-3 


0.3661e-3 


0.3601e-3 


2.15 


0.3583e-3 




^12 


0.3927e-3 


0.3661e-3 


0.3601e-3 


2.15 


0.3583e-3 




OJ22 


0.5983e-3 


0.5446e-3 


0.5321e-3 


2.11 


0.5284e-3 




0)11 


0.1772e-3 


0.1760e-3 


0.1757e-3 


2.10 


0.1756e-3 




0)21 


0.3634e-3 


0.3592e-3 


0.3584e-3 


2.33 


0.3582e-3 




0)12 


0.3650e-3 


0.3594e-3 


0.3584e-3 


2.53 


0.3582e-3 




o>22 


0.5440e-3 


0.5312e-3 


0.5286e-3 


2.30 


0.5279e-3 




0)11 


0.1779e-3 


0.1761e-3 


0.1757e-3 


1.94 


0.1755e-3 




0)21 


0.3654e-3 


0.3599e-3 


0.3585e-3 


1.95 


0.3580e-3 




^12 


0.3679e-3 


0.3600e-3 


0.3585e-3 


2.40 


0.3582e-3 




OJ22 


0.5489e-3 


0.5325e-3 


0.5289e-3 


2.18 


0.5278e-3 



It can be seen from Table 5 that as for the source problem, our Method 2 lead 
to wrong result for triangular meshes when the thickness of the plate is small, see 
Remark 3. For any other family of meshes the method is locking free and converges 
with a quadratic order. 



Table 6 shows the four lowest vibration frequencies computed by Method 2 with 
successively refined meshes of each type for a simply supported plate with thickness t = 
0.01. The table includes orders of convegence, as well as accurate values extrapolated 
by means of a least-squares fitting. Furthermore, the last two columns show the results 
reported in [30, 25]. In every case, we have used a Poisson ratio v = 0.3 and a correction 
factor k = 0.8333. The reported non-dimensional frequencies are independent of the 
remaining geometrical and physical parameters, except for the thickness-to-span ratio. 
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Table 6 Lowest non-dimensional vibration frequencies for a SSSS square plate and t = 0.01 



Mesh 


Mode 


N = 16 


N = 32 


N = 64 




TP/vf ran 

J_j Jv L 1 CLJJ . 


ho 


25 




LOW 


0966 


0963 

V./ . V7 -7 V ' '7 


0963 


2.04 


0963 


0963 

\ ) . V7 .7 \J it 


0963 


T- 


L021 


2416 


U.Z4UO 


n 9 mfi 




O 940R 
U. Z^iUU 


940R 


9A0fi 




LO\2 


0.2425 


0.2411 


0.2407 


2.02 


0.2406 


0.2406 


0.2406 




L022 


u.oooy 


U.OOUO 


U.OOJU 


9 nn 


O "\9.A7 
U.OO'l 1 


^8.47 


U.OOUO 




LOW 


0.0967 


0.0964 


0963 


1.96 


0963 


0963 

V7 . V7 -7 17 'J 


0963 

V7 . V7 -7 V7 ' 7 


>h 


Ld21 


0.242 4 


0.2410 


0.2407 


1.91 


0.2406 


0.2406 


0.2406 




LO\2 


0.2434 


0.2413 


0.2408 


1.93 


0.2406 


0.2406 


0.2406 




L022 


0.3914 


n 3Kfi^ 

U.OOUO 


n 38^9 


1.92 


0.3847 


0.3847 


0.3848 




LOW 


0966 


0.0964 


0963 


2.00 


0963 


0963 

V ' < V7 - 7 V7 ' J 


0963 

v 7 . y.t >y \J ' t 


T 4 
>h 


L02Y 


2426 


2411 


fl 9407 


2 02 




n 9zLnfi 




LOV1 


2426 


2411 


fl 9407 


2 02 










L022 


3SQ8 
u . ooyo 


u. oouu 


U. OOUU 


2.01 


0.3847 


0.3847 


0.3848 




LOW 


0.0967 


0.0964 


0.0963 


2.01 


0.0963 


0.0963 


0.0963 


■T~> 
' h 


L02\ 


0.2429 


0.2411 


0.2407 


2.07 


0.2406 


0.2406 


0.2406 




LO\2 


0.2441 


0.2414 


0.2408 


2.09 


0.2406 


0.2406 


0.2406 




L022 


0.3910 


0.3863 


0.3851 


2.01 


0.3847 


0.3847 


0.3848 




LOW 


0.0964 


0.0963 


0.0963 


2.82 


0.0963 


0.0963 


0.0963 


7? 

' h 


L02\ 


0.2411 


0.2406 


0.2406 


3.79 


0.2406 


0.2406 


0.2406 




LO\2 


0.2417 


0.2407 


0.2406 


3.82 


0.2406 


0.2406 


0.2406 




L022 


0.3869 


0.3850 


0.3848 


3.40 


0.3848 


0.3847 


0.3848 




LOW 


0.0965 


0.0963 


0.0963 


2.53 


0.0963 


0.0963 


0.0963 


T~ 
>h 


L02\ 


0.2416 


0.2408 


0.2406 


2.35 


0.2406 


0.2406 


0.2406 




LO\2 


0.2427 


0.2408 


0.2407 


3.48 


0.2406 


0.2406 


0.2406 




L022 


0.3889 


0.3854 


0.3848 


2.61 


0.3847 


0.3847 


0.3848 



Table 7 shows the four lowest vibration frequencies computed by Method 2 with 
successively refined meshes of each type for a plate with a free edge (with three clamped 
edges and the fourth free) with thickness t = 0.01. The table includes orders of con- 
vegence, as well as accurate values extrapolated by means of a least-squares fitting. 
Furthermore, the last two columns show the results reported in [30,25]. In every case, 
we have used a Poisson ratio v = 0.3 and a correction factor k = 0.8601. The reported 
non-dimensional frequencies are independent of the remaining geometrical and physical 
parameters, except for the thickness-to-span ratio. 
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Table 7 Lowest non-dimensional vibration frequencies for a CCCF square plate and t = 0.01 



Mesh 


Mode 


— 32 


N = 64 


TV = 128 




TP/vf ran 

J_J^V L 1 CLJJ ■ 


ho 

[OU] 


25 

[ZD] 






0.1215 


0.1179 


0.1169 


1.98 


0.1166 


0.1166 


0.1171 


T 4 
>h 




u. zuou 


n i Q7n 


u. lyji 


1 Q7 


n 1 QAQ 
u. lyyy 


n 1 QAQ 


u. lyji 




^12 


u.oooo 


0.3144 


u.ouyu 


2.14 


0.3081 


u. ouoo 


u.ouyo 




U22 




U.J / IO 


u.o / 


1 Q7 


u.o / ou 


u.o t ou 


u.o / ^±u 




UJ\ 1 


0.1199 


0.1173 


0.1167 


2.18 


0.1166 


0.1166 


0.1171 


n 


W21 


0.1986 


0.1958 


0.1951 


2.00 


0.1948 


0.1949 


0.1951 




wia 


0.3264 


0.3117 


0.3086 


2.25 


0.3078 


0.3083 


0.3093 




^22 


0.3791 


0.3749 


0.3738 


1.92 


0.3734 


0.3736 


0.3740 




wn 


0.1177 


0.1169 


0.1167 


2.00 


0.1166 


0.1166 


0.1171 




CJ 21 


0.1967 


0.1953 


0.1949 


2.13 


0.1948 


0.1949 


0.1951 


U>12 


0.3134 


0.3090 


0.3081 


2.22 


0.3079 


0.3083 


0.3093 




1^22 


0.3753 


0.3738 


0.3735 


2.30 


0.3734 


0.3736 


0.3740 




wn 


0.1180 


0.1169 


0.1167 


1.90 


0.1166 


0.1166 


0.1171 


n 


^21 


0.1974 


0.1954 


0.1950 


2.15 


0.1948 


0.1949 


0.1951 




wia 


0.3151 


0.3095 


0.3082 


2.08 


0.3078 


0.3083 


0.3093 




1^22 


0.3772 


0.3743 


0.3736 


2.22 


0.3734 


0.3736 


0.3740 



It can be seen from Tables 6 and 7 that our method converges with a quadratic 
order. 



4.3 Buckling of plates 

The effectiveness of the MDF method for buckling analysis are demostrated by exam- 
ples with different thickness, boundary conditions and different in-plane compressive 
stress cr. 

We have computed approximations of the buckling coefficients X bc = \ bp t 2 being 
the smallest (the critical load) by which the chosen in-plane compressive stress cr must 
be multiplied by in order to cause buckling. In order to compare our results with those 
in [33,34,42], a non-dimensional buckling intensity is defined as: 

\ bc T 

7T 2 D' 

here \ b h c = \ hp t 2 are the computed buckling coefficients, L is the plate side length and 
D is the flexural rigidity defined as D = £i 3 /[12(l - v 2 )}. 

4-3.1 Uniformly compressed plate 

In this couple of tests, we use cr = I, corresponding to a uniformly compressed plate 
(in the x, y directions). 

First, we consider a simply supported plate, since analytical solutions are available 
(see [45]) for that case. In Table 8, we report the four lowest non-dimensional buckling 
intensities K\, . . . , K4, for the thickness t = 0.01, and L = 1 computed by Method 3 
with four different family of meshes. The table includes computed orders of convergence, 
as well as more accurate values extrapolated by means of a least-squares procedure. 
Furthermore, the last column reports the exact buckling intensities. In this case, we 
have used a Poisson ratio v = 0.3 and a correction factor k = 5/6. 
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Table 8 Lowest non-dimensional buckling intensities K\ , . . . , K4 for a SSSS square plate and 
t = 0.01 



Mesh 


K 


N = 16 


N = 32 


N = 64 


Order 


Extrap. 


Exact 




Ki 


2.0315 


2.0075 


2.0011 


1.90 


1.9987 


1.9989 


' h 


1<2 


5.1607 


5.0381 


5.0046 


1.87 


4.9920 


4.9930 




K 3 


5.2262 


5.0541 


5.0086 


1.92 


4.9922 


4.9930 




K 4 


8.5170 


8.1249 


8.0186 


1.88 


7.9788 


7.9820 






2.0381 


2.0086 


2.0013 


2.01 


1.9989 


1.9989 


r 4 

' h 


1<2 


5.2116 


5.0465 


5.0063 


2.04 


4.9934 


4.9930 




K 3 


5.2116 


5.0465 


5.0063 


2.04 


4.9934 


4.9930 




K 4 


8.6292 


8.1388 


8.0209 


2.06 


7.9839 


7.9820 




K x 


2.0412 


2.0093 


2.0015 


2.02 


1.9989 


1.9989 


'h 


K 2 


5.2242 


5.0485 


5.0063 


2.06 


4.9931 


4.9930 




K 3 


5.2929 


5.0641 


5.0099 


2.08 


4.9932 


4.9930 




K 4 


8.6788 


8.1504 


8.0234 


2.06 


7.9836 


7.9820 




K x 


2.0347 


2.0068 


2.0010 


2.27 


1.9995 


1.9989 


'h 


K 2 


5.1962 


5.0424 


5.0053 


2.05 


4.9935 


4.9930 




K 3 


5.2429 


5.0451 


5.0063 


2.35 


4.9968 


4.9930 




K 4 


8.5851 


8.1192 


8.0158 


2.17 


7.9862 


7.9820 



It can be seen from Table 8 that our method converges to the exact values with a 
quadratic order. 



As a second test, we present the results for the lowest non-dimensional buckling 
intensity K\ for a clamped plate with varying thickness t, in order to assess the stability 
of the Method 3 when t goes to zero. It is well known that K\ converges to the non- 
dimensional buckling intensity of an identical Kirchhoff-Love uniformly compressed 
clamped plate. 



In Table 9, we report the lowest non-dimensional buckling intensity K\ of a uni- 
formly compressed clamped plate with varying thickness t and L = 1. We have used 
five different family of meshes. The table includes computed orders of convergence, as 
well as more accurate values extrapolated by means of a least-squares procedure. In the 
last row of each family of meshes we report the limit values as t goes to zero obtained 
by extrapolation. In this case, we have used a Poisson ratio v = 0.25 and a correction 
factor k = 5/6. 
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Table 9 Lowest non-dimensional buckling intensity K\ of a clamped plate with varying thick- 
ness. 



Mesh 


t 


N = 16 


N = 32 


N = 64 


Order 


Extrap. 




0.1 


4.9031 


4.6658 


4.6099 


2.09 


4.5929 


'h 


0.01 


7.1314 


5.5307 


5.3275 


2.98 


5.2982 


0.001 


9.9817e+l 


9.2546 


5.6315 


4.00 


4.3035 




0.0001 


9.2723e+3 


2.6878e+2 


1.2923e+l 


4.00 


-0.1678 




(extrap.) 














0.1 


4.6316 


4.6015 


4.5937 


1.93 


4.5909 




0.01 


5.3423 


5.3072 


5.2981 


1.96 


5.2950 


0.001 


5.3509 


5.3157 


5.3066 


1.96 


5.3035 




0.0001 


5.3510 


5.3158 


5.3067 


1.96 


5.3035 




(extrap.) 


5.3510 


5.3158 


5.3067 


1.96 


5.3036 




0.1 


4.6476 


4.6058 


4.5948 


1.92 


4.5908 




0.01 


5.3611 


5.3121 


5.2994 


1.94 


5.2949 




0.001 


5.3697 


5.3207 


5.3079 


1.94 


5.3034 




0.0001 


5.3698 


5.3208 


5.3080 


1.94 


5.3035 




(extrap.) 


5.3698 


5.3208 


5.3080 


1.94 


5.3035 




0.1 


4.6441 


4.6043 


4.5943 


2.00 


4.5910 


r 4 


0.01 


5.3564 


5.3103 


5.2989 


2.01 


5.2951 




0.001 


5.3649 


5.3188 


5.3074 


2.01 


5.3036 




0.0001 


5.3649 


5.3189 


5.3074 


2.01 


5.3037 




(extrap.) 


5.3650 


5.3189 


5.3074 


2.01 


5.3037 




0.1 


4.6487 


4.6054 


4.5946 


2.00 


4.5909 


r 5 

'h 


0.01 


5.3702 


5.3128 


5.2993 


2.09 


5.2952 




0.001 


5.3863 


5.3240 


5.3085 


2.01 


5.3034 




0.0001 


5.3866 


5.3242 


5.3088 


2.01 


5.3036 




(extrap.) 


5.3867 


5.3242 


5.3087 


2.01 


5.3036 



Additionally, we have also computed the lowest buckling intensity of a Kirchhoff- 
Love plate by using the finite element method analyzed in [41]. 

In Table 10, we report the lowest non-dimensional buckling intensity of a uniformly 
compressed clamped plate with L = 1. In this case we considered a Poisson ratio 
v = 0.25. 



Table 10 Lowest non-dimensional buckling intensity of a uniformly compressed clamped thin 
plate (Kirchhoff-Love model) computed with the method from [41]. 



Method 


N = 24 N = 36 N = 48 N = 60 Order 


Extrapolated 


[41] 


5.3051 5.3042 5.3039 5.3038 2.61 


5.3037 



It is clear from the results of Tables 9 and 10, that our Method 3 lead to wrong 
result for triangular meshes when the thickness of the plate is small, see Remark 3. 
For all the other family of meshes the method is locking free and do not deteriorate as 
the plate thickness become smaller. 
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4-3.2 Plate uniformly compressed in one direction 



In this couple of tests, we use 



1 




corresponding to a plate subjected to uniaxial compression (in the x direction). We 
consider different boundary conditions. 

In Tables 11 and 12, we report the lowest non-dimensional buckling intensity K\, for 
a clamped and simply supported plate, respectively, with thickness t = 0.1, and L = 1 
computed by Method 3 with different family of meshes. The table includes computed 
orders of convergence, as well as more accurate values extrapolated by means of a 
least-squares procedure. Furthermore, the last two columns show the results reported 
in [33,34]. In these cases, we have used a Poisson ratio v = 0.3 and a correction factor 
k = 5/6. 



Table 11 Lowest non-dimensional buckling intensity K\ for a CCCC square plate and t = 0.1 



Mesh 


K 


N = 32 


N = 64 


N = 128 


Order 


Extrap. 


[33] 


N 


T' A 


Ki 


8.3849 


8.3157 


8.2978 


1.95 


8.2915 


8.2917 


8.2931 


j% 

'h 


Ki 


8.4222 


8.3260 


8.3004 


1.91 


8.2912 


8.2917 


8.2931 


r 4 


Ki 


8.3987 


8.3185 


8.2984 


2.00 


8.2917 


8.2917 


8.2931 




K 1 


8.4273 


8.3255 


8.3001 


2.00 


8.2916 


8.2917 


8.2931 


'h 


Ki 


8.3663 


8.3110 


8.2963 


1.91 


8.2910 


8.2917 


8.2931 






8.3715 


8.3121 


8.2965 


1.93 


8.2909 


8.2917 


8.2931 



Table 12 Lowest non-dimensional buckling intensity Ki for a SSSS square plate and t = 0.1 



Mesh 


K 


N = 32 


N = 64 


N = 128 


Order 


Extrap. 


[33] 


[34] 


T' 2 


Ki 


3.7993 


3.7897 


3.7873 


1.97 


3.7864 


3.7865 


3.7873 


T' 1 


Ki 


3.8029 


3.7907 


3.7875 


1.96 


3.7864 


3.7865 


3.7873 




K 1 


3.8049 


3.7911 


3.7876 


2.00 


3.7864 


3.7865 


3.7873 




Ki 


3.8062 


3.7913 


3.7877 


2.01 


3.7865 


3.7865 


3.7873 




K 1 


3.7975 


3.7892 


3.7871 


1.98 


3.7864 


3.7865 


3.7873 


>h 


K x 


3.8011 


3.7903 


3.7874 


1.90 


3.7863 


3.7865 


3.7873 



It can be seen from Tables 11 and 12 that our method converges with a quadratic 
order. 

4-3.3 Shear loaded plate 
In this test, we use 

"0 1" 
a ~ [10_ ' 

corresponding to a plate subjected to shear load. We consider different boundary con- 
ditions. 
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In Table 13, we report the lowest non-dimensional buckling intensity K\. for a 
simply supported plate with thickness t = 0.01, and L = 1 computed by Method 3 
with different family of meshes. The table includes computed orders of convergence, 
as well as more accurate values extrapolated by means of a least-squares procedure. 
Furthermore, the last two columns show the results reported in [42]. In these cases, we 
have used a Poisson ratio v = 0.3 and a correction factor k = 5/6. 

Table 13 Lowest non-dimensional buckling intensity K\ for a SSSS square plate and t = 0.01 



Mesh 


K 


TV = 32 


N = 64 


N = 128 


Order 


Extrap. 


[42] 


T' 2 


Ki 


9.4832 


9.3514 


9.3180 


1.98 


9.3067 


9.2830 


T*- 


Kr 


9.3848 


9.3270 


9.3119 


1.94 


9.3066 


9.2830 


Tt- 


Ki 


9.4602 


9.3450 


9.3164 


2.01 


9.3069 


9.2830 




Kj 


9.4759 


9.3483 


9.3170 


2.03 


9.3069 


9.2830 


r^- 


K X 


9.4196 


9.3346 


9.3134 


2.01 


9.3064 


9.2830 


'h 


A'i 


9.4640 


9.3460 


9.3163 


1.99 


9.3063 


9.2830 



It can be seen from Table 13 that our method converges with a quadratic order. 



5 Conclusions 

We assessed numerically the actual performance of the method proposed in [11], extend- 
ing it also to free vibration and buckling problems of plates. We tested different families 
of mimetic meshes, different values of the relative thickness and various boundary con- 
ditions. In all the three types of problems considered (source problem, free vibration, 
buckling) the method was shown to be locking free and to converge with an optimal 
rate both in discrete L°° and H 1 norms for meshes made with elements with 4 or more 
edges. In some occasions, a super convergence rate was noticed. Moreover, differently 
from standard quadrilateral finite elements, the method shows a robust behavior also 
for uniformly distorted families of meshes such as those in Figure 5. We thus conclude 
that the proposed method is very reliable for Reissner-Mindlin plate computations. 

Acknowledgements The third author was partially supported by CONICYT-Chile through 
FONDECYT project No. 11100180, and by Centro de Investigation en Ingenierfa Matematica 
(CI 2 MA), Universidad de Conception, Chile. 
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